function V = feature_vquan(imgs, feafun, K0, K1, p)
%FEATURE_VQUAN Feature vector quantization
%
%   V = FEATURE_VQUAN(imgs, feafun, K0, K1);
%   V = FEATURE_VQUAN(imgs, feafun, K0, K1, p);
%
%       Use two-level K-means to quantize feature vectors from images,
%       returing a collection of center vectors.
%
%       The procedure consists of two stages:
%       - clustering the feature vectors from each image, getting K0
%         centers.
%       - The center vectors from all images are pooled, which are then
%         clustered into K1 centers.
%
%       Input arguments:
%       - imgs:     The cell array of images
%
%       - feafun:   The function that transforms an image into a matrix
%                   comprised of feature vectors.
%
%                       X = feafun(I);
%
%                   This statement should return a matrix of size d x n,
%                   where d is the feature dimension, and n is the number
%                   of feature vectors.
%
%       - K0:       The number of centers for each image.
%
%       - K1:       The size of final dictionary.
%                   pre-condition: K1 < K0 * numel(imgs).
%
%       - p:        sub-sampling ratio. (p <= 1)
%
%
%       Output arguments:
%       - V:        The resultant center vectors, size: d x K1.
%

%% verify input arguments

if ~iscell(imgs)
    error('feature_vquan:invalidarg', 'imgs should be a cell array.');
end
n = numel(imgs);

if ~isa(feafun, 'function_handle')
    error('feature_vquan:invalidarg', 'feafun should be a function handle.');
end

if ~(isnumeric(K0) && isscalar(K0) && K0 == fix(K0) && K0 >= 1)
    error('feature_vquan:invalidarg', 'K0 should be a positive integer scalar.');
end

if ~(isnumeric(K1) && isscalar(K1) && K1 == fix(K1) && K1 >= 1)
    error('feature_vquan:invalidarg', 'K1 should be a positive integer scalar.');
end

if K1 >= K0 * n
    error('feature_vquan:invalidarg', ...
        'K1 should have K1 < K0 * numel(imgs).');
end

if nargin < 5
    p = 1;
end


%% main

disp('Stage 1: Per-image clustering ...');

Cs = cell(1, n);
ws = cell(1, n);

for i = 1 : n    
    fprintf('\ton image %d / %d\n', i, numel(imgs));
    X = feafun(imgs{i});    
    
    if p < 1
        nx = size(X, 2);
        X = X(:, rand(1, nx) < p);
    end
    
    [L, Cs{i}] = kmeans_std(X, K0, ...
        'MaxIter', 30, 'Display', 'off', 'Init', 'random');      
    w = intcount(K0, L);
    ws{i} = w / sum(w);
end

disp('Stage 2: overall clustering ...');

Cs = [Cs{:}];
ws = [ws{:}];

[~, V] = kmeans_std({Cs, ws}, K1, ...
    'MaxIter', 300, 'Display', 'off', 'Init', 'km++');


